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We report an order-A approach to compute the Kubo Hall conductivity for disorderd two- 
dimensional systems reaching tens of millions of orbitals, and realistic values of the applied external 
magnetic fields (as low as a few Tesla). A time-evolution scheme is employed to evaluate the Hall 
conductivity a xy using a wavepacket propagation method and a continued fraction expansion for the 
computation of diagonal and off-diagonal matrix elements of the Green functions. The validity of 
the method is demonstrated by comparison of results with brute-force diagonalization of the Kubo 
formula, using (disordered) graphene as system of study. This approach to mesoscopic system sizes 
is opening an unprecedented perspective for so-called reverse engineering in which the available 
experimental transport data are used to get a deeper understanding of the microscopic structure 
of the samples. Besides, this will not only allow addressing subtle issues in terms of resistance 
standardization of large scale materials (such as wafer scale polycrystalline graphene), but will also 
enable the discovery of new quantum transport phenomena in complex two-dimensional materials, 
out of reach with classical methods. 

PACS numbers: 73.43.-f,72.80.Vp, 73.22.Pr 


I. INTRODUCTION 

The Hall effect is one of the central phenomenon in 
Condensed Matter physics with great importance for 
measuring carrier density and mobility in semiconduc¬ 
tors. The long history started almost 150 years ago when 
E. H. Hall observed a transverse voltage drop over a con¬ 
ducting bar driven by a longitudinal current^ and it expe¬ 
rienced an important turn in 1980 when K. von Klitzing 
observed a quantized version of this effect, the so-called 
Quantum Hall Effect (QHE)4 This had an immediate 
impact and triggered huge amount of work which led to 
the definition of a resistance standard^ 

The effect of disorder and related localization phe¬ 
nomenon is central to understand the integer QHE from a 
bulk perspective, where a xy plateaus develop by varying 
the charge density, while the plateau region coincide with 
a region of vanishing o xx 4 This is explained by Anderson 
localization of the wave functions in the Landau levels in¬ 
duced by disorder. The most standard method to treat 
bulk conductivities microscopically is the linear-response 
theory, or the Kubo formula, which was first applied by 
Aoki and Andc4 to the QHE problem^. Such quantity 
can be connected to measurements in the Corbino geom¬ 
etry for which electrodes are attached to the inner and 
outer perimeters of an annular sample, which allow to 
investigate bulk transport physics in the high magnetic 
field regime. 

Recently, the QHE has played a crucial role for the 
first unambigous proof of the existence of single atomic 
layers of grapheneSrs. Indeed, in graphene, the pecu¬ 
liar Berry phase related to the pseudospin quantum de¬ 
gree of freedom results in a unique spectrum with four¬ 
fold degenerate Landau levels (due to the spin and valley 


degeneracies), and half-integer QHE with Hall plateaus 
emerging at a xy = ^-(n + 1/2), with integers n. QHE 
related research in graphene is a vivid field, study¬ 
ing electron-electron interactional^!!, superlattice44r— , 
graphene functionalization and topological currents. Ad¬ 
ditional transverse transport phenomena are being stud¬ 
ied such as (quantum-) anomaleous Hall effect ) 15 : 16 , spin 
Hall effect— and quantum spin Hall effect— 

The role of disorder in graphene is also debated 
and the QHE features, with very robust behavior to 
large disorder—, also present additional peculiarities de¬ 
pending on the symmetry breaking aspects conveyed 
by defects— — Recently, we have studied the case of 
oxidized graphene, and have shown the formation of 
disorder-induced resonant critical states appearing in the 
zero-energy Landau level with finite a xx and suggesting 
a zero-energy a xy quantized plateau—. The confirmation 
of such a plateau will demand to revise the description of 
topological invariant in disordered graphene, an exciting 
direction of work. However a real space calculation of 
a xy is needed and demand for the development of new 
methodology. 

Indeed, the importance of an efficient computational 
methodology for the Hall conductivity a xy stands in 
sharp contrast to the limited possibilities to simulate this 
quantity with usual approaches. One of them starts with 
the Kubo formula and requires the knowledge of the full 
eigenstates of a given system^. Such formulation allows 
for advanced analysi s 26 ’ 27 but suffers from unfavourable 
cubic scaling with system size which restricts to small 
systems because diagonalization is necessary. All eigen¬ 
states are needed and have to be combined to calculate 
a xy . This computational limitation makes the analysis 
of realistic models of disordered samples a priori im- 
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possible. Additionally, such approach is restricted to 
unreasonably strong magnetic fields that exceed by far 
what is achievable experimentally, and forces the mag¬ 
netic length scale to dominate over all other length scales 
of the problem (mean free path, average distance between 
impurities,...) limiting the exploration of complex mag¬ 
netotransport in intermediate or even low-magnetic field 
regimes, thus calling for approaches with improved nu¬ 
merical performance—r— 

In this paper, we present a highly efficient order-N 
algorithm for the computation of the Kubo Hall con¬ 
ductivity that allows us to circumvent aforementioned 
limitations, and offer fascinating perspectives for explor¬ 
ing transverse transport phenomena in disordered two- 
dimensional materials with almost arbitrary complexity. 
A validation of the method on several models of clean or 
disordered graphene is made by comparing this new time 
evolution method with brute force cliagonalization tech¬ 
nique. Some preliminary illustration of the method have 
been published in a prior work— The method is inspired 
by the real space implementation of the dissipative Hall 
conductivity at zero frequency, which has proven undis¬ 
puted computational efficiency and predictive power (see 
Ref i 31 i 32 for computational details, and Ref— for some 
illustration on realistic models of disordered graphene). 


II. METHODOLOGY 

Kubo formalism - The general framework is provided 
by Kubo’s linear response theor y 34 ' 35 in which the dc 
conductivity is given as 

l r°° rP 

<j xy {u = 0 ) = — / dt dXTi [ pojyj x (t + ih A)] ( 1 ) 
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has been widely employed in the literature for studying 
model systems with phenomenological description of dis¬ 
order. Unfortunately such numerical operations quickly 
makes the methodology computationally prohibitive, in 
particular for simulating complex and large system sizes 
closer to the experimental ones and at moderate mag¬ 
netic fields. Clearly, the computational problem then 
becomes untreatable with diagonalization-based meth¬ 
ods, i.e. it becomes impossible to calculate (and store) 
all eigenvalues and eigenvectors, so that time-evolution 
based approaches appear promising. By introducing 
j_ oo dE8{E — H) we arrive at an alternative represen¬ 
tation of the Hall conductivity 

2 AC r°° r°° 

<7xy =- y J dt J dEf(E - C) 
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, (4) 

which proves useful for the numerical implementation. 
In Eq. the trace is replaced by an average over 

a random-phase state |'Ri) according to JO fc (fc|...|fc) —t 
.ZV s (\f r i|...|\I>i). Such states are normalized by using the 
number of sites N s . We further introduce the projection 
operator J2f=i as well as the time-dependent 

cr(t) fulfilling a xy = lim t _> 0 0 a xy (t) with 


lim Re 

r)—, >0+ 


<J X y{t ) 


2 

7 lV s 



/ OO 

dEf(E - C) 

-OO 


Nr 

lim Vim [«,(£)] 
7}->o+ J 

j= 1 


Re 


(*j| 3yU\t) \ . j x U{t)\*i) 

th — ri y IT] 


( 5 ) 


Starting from Eq. 0. the evaluation of the conductiv¬ 
ity &xy proceeds by introducing an orthonormal many- 
particle basis with (n\H\n) = E n which allows to rewrite 
Eq. 0 into 


Oxy 
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where i] is a small parameter necessary to converge this 
expression. 

In the case of non-interacting electrons one can proceed 
further towards the single-particle picture and obtain 
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where the Fermi-Dirac distribution is described by /(fffc — 
C) and the chemical potential £ has been introduced. 
The single-particle eigenvalues £k can in principle be ob¬ 
tained by matrix diagonalization of finite systems which 


where 


*A E) = (6) 

are the the matrix elements of the Green function (with 
z = E + ir]) and V s = V/N s is the volume per site. 
Equation ® is the basis for the numerical evaluation of 
the Hall conductivity for which one eventually takes the 
limit tf —> oo. One notes, however, that the result does 
not converge in all cases when taking this limit and an os¬ 
cillatory behaviour of the integral is generally observed, 
particularly for the case of clean systems in high mag¬ 
netic fields (and well separate Landau levels), which is 
the usual test reference for Kubo Hall transport simula¬ 
tions. 

The reason is that the Kubo formula Eq. 0 assumes 
that at “infinite past” times, the system is in equilib¬ 
rium, i.e. at zero electric field, the current should vanish. 
In other words, the current-current correlation function 
should decay with time difference t. This should be ful¬ 
filled also for the final result. In pristine systems this is 
not the case and must be imposed a posteriori. In order 













to take this into account an additional small parameter 
S > 0 has to be introduced such that the integral 
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converges. 

We note that the same parameter has to be included in 
full analogy in the standard expression for the calculation 
of the Hall conductivity based upon the knowledge of the 
eigensystem. We here display such formula^^ 
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which is identical to Ref^ for the reasonable choice that 

= S. These two equations 0 and 0 are comparable 
and will be evaluated for a couple of examples in Sect, 
mu as numerical checks. 

Numerical implementation - The calculation proceeds 
by splitting the time-dependent and time independent 
parts ( Kj(E )) while the sum over j in 0 controls the en¬ 
ergy dependent convergence (see details below). Firstly, 
in order to evaluate Eq. we make use of the Lanc- 
zos tridiagonalization of the Hamiltonian H which allows 
us to use the continued fraction expansion for the diag¬ 
onal Green function. From the diagonal elements one 
can obtain the necessary off-diagonal elements in a sec¬ 
ond iterative step. We find a simple recursion relation 
connecting both which reads for n > 1 

Kn-\-1 — j [ bn — lKn—1 Ura)^n] (9) 

On 

while the first step in the iteration is given as K 2 = 
FT [1 - 

Secondly, the time evolutions of the quantities | Ti} and 
Ijy'&j) can be performed efficiently through a Cliebychev 
expansion of the time-evolution operator U (t) which usu¬ 
ally converges quite rapidly. For the final calculation of 
the brackets in Eq. 0 we have to calculate off diago¬ 
nal matrix elements with the advanced Green operator 
(z — H)- 1 . For the calculation of this second set of off- 
diagonal Green functions in Eq. 0, we use the continued 
fraction expansion. All operations are linear in sample 
size and this holds also for the total expression since the 
scalar product in j is restricted to a cutoff value Nr. 
We use S = 0.002 if not stated otherwise throughout the 
paper. 

System description The electronic properties of 
graphene can be described by the standard orthogonal 



FIG. 1: (color online) Absolute values of the off-diagonal ma¬ 
trix elements Eq. @ for pristine graphene (a) and disordered 
graphene (b). Different curves correspond to energies around 
the lowest-energy Landau level. 


tight-binding model 

U = Y j y a W){a\- 70 £ e-^|a)(/3|, (10) 

a ( a ,P) 

where 70 = 2.7eU is the nearest neighbor transfer inte¬ 
gral and V a the on-site energy, which can be chosen to 
describe various disorder models. Spatially uncorrelated 
Anderson disorder is introduced through on-site energies 
taken at random from [—W 70 / 2 , W 70 / 2 ] where W gives 
the disorder strength. This is a commonly used disor¬ 
der model for exploring the metal-insulator transition 
in low-dimensional system s 26 ’ 37 , which, in average, does 
not break electron-hole symmetry or inversion symmetry. 
The inversion symmetry of the system related to the sym¬ 
metry of AB sublattices can be artificially broken by a 
staggered potential of the form V a = +(—)Vab for A(B) 
sublattices (including all atoms). A less artificial vari¬ 
ant of such potential can be introduced by diluting such 
sublattice terms for a weaker and random distribution, 
which mimicks some correlation in the potential beyond 
the Anderson model. 

The magnetic field is implemented through a Peierls 
phased added to 70 which determines the magnetic flux 
to <l) = § A ■ dl = h/e Ehexagon P er plaquette.^ 
Spin degeneracy is assumed throughout the paper and 
an anti-symmetrization procedure for a xy (E) has been 
performed consistent with electron-hole symmetry for all 
considered disorder potentials. 


III. RESULTS 

Before we enter into the discussion of the physical re¬ 
sults, we illustrate the convergence behaviour of the nu¬ 
merical algorithm which can be controled by the expres¬ 
sion ©. A key quantity for the convergence is Kj{E) 
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FIG. 2: (color online) Hall conductivity for varying magnetic 
field strength. Main frame: moderately high fields; lower 
inset: extremely fields. Dashed curves in all frames indicate 
simulations with time-evolution Kubo (TEK) approach. Note 
the magnetic-field scaling of Landau level energies E oc y/B 
which is reflected in varying energy scales of main frame and 
inset (70 units). 


FIG. 3: (color online). Hall conductivity for Anderson disor¬ 
der (dashed curves: exact diagonalization; solid curves: TEK) 
at different magnetic fields (main frame). Effect of ncreasing 
disorder for high field 963T (upper inset) and intermediate 
fields of 45T (lower inset). 


which can be studied prior to the time evolution of 
wavepackets and allows to estimate the required number 
of recursion steps which depend on the specific physical 
details (e.g. disorder, magnetic field), on the selected 
energy, and on the broadening 77. Fig. |T) presents the 
modulus of k plotted versus j for some selected ener¬ 
gies in the case of pristine graphene (a) and disordered 
graphene (b). As a general trend, Fig. |T] shows that in 
case of pristine systems, the convergence of the Lanczos 
recursion can be reached with a low number of recursion 
steps (Nr), while increasing disorder requires Nr to be 
increased typically to a few thousands. 

Turning to the simulations of the Hall conductivity, 
we first discuss the results for pristine graphene in Fig. 
[ 2 ] which shows the comparison of diagonalization method 
and time-evolution Kubo approach. Figure 1 (inset) 
shows the case of extremely high magnetic fields (fluxes). 
We recall that, in order to fulfill the boundary conditions 
of the periodic system, the limit of very high fields is usu¬ 
ally considered by diagonalization-based studies, because 
only small sample sizes can be treated by matrix diag¬ 
onalization. The magnetic field is simultaneously given 
in the inset of Fig. [2] in terms of the total flux pene¬ 
trating the sample. The corresponding results from the 
time-evolution Kubo method (TEK) is plotted as dashed 
lines for the same magnetic fields (values indicated in 
the inset legend). For the selected energies, the curves 
coincide visually. The square root scaling of their posi¬ 
tion with magnetic field is evident and Hall-conductance 
steps occur at the expected positions and with the ex¬ 
pected height corresponding to the half-integer sequence 
a xy = ±( 1/2 + n)4e 2 /h. 

Next, we consider the interesting case of moderately 


high magnetic fields such as frequently observed in cur¬ 
rent experiments^— in the main frame of Fig. [2] The 
Kubo approach allows to reduce magnetic fields down to 
only few Teslas where the first Landau levels appear at 
Fermi energies of few tens of meV. Still the plateau se¬ 
quence is obtained with very good accuracy, which shows 
that such approach can be of great practical use. Most 
notably, the figure demonstrates conductance quantiza¬ 
tion for fields as low as 11 Tesla. As an effect of finite 
energy resolution, we observe that the relative sharpness 
of the steps is reduced with lowering the field in this 
regime. The conductance slope starts to become noti- 
cable as Landau levels get closer (as seen for 11.4 T at 
<j xy = ±10 e 2 /h). This consequence of the considered en¬ 
ergy resolution (here 77 = 0.00025) can be systematically 
improved further by increasing the number of recursion 
steps and decreasing 77 . 

As a second example to illustrate the performance of 
the algorithm, we focus on the case of homogeneously 
disordered graphene. The case of Anderson disorder is 
chosen for Fig. [3] because this disorder is expected to 
induce a transition from the QHE regime to the con¬ 
ventional Anderson insulating state for sufficiently high 
value of W. The main-frame indicates very good agree¬ 
ment (within few percents) obtained between the Kubo 
algorithm and exact diagonalization techniques at low 
energy for very high magnetic fields. At higher energy, 77 
allows for fast convergence, at the expense of small loss in 
information compared to the exact method when Landau 
levels are getting closer in energy. To illustrate the effect 
of increasing disorder, a zoom on the first two Hall steps 
is provided for 963 T in the upper inset. The Landau lev¬ 
els are broadened with increasing W around the critical 
states, as expected from perturbation theory. Simulta¬ 
neously, the first and second plateaus remain at ± 2 e 2 /h 
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FIG. 4: (color online) Hall conductivity for p=2.5% of AB- 
sublattice breaking defects (strenght V s = O. 270 ) distributed 
at random in space (see main text). Vertical dotted line indi¬ 
cates nominal gap of a correspondingly homogenous AB po¬ 
tential of strength pVab = O.OO 570 . 


and ± 6 e 1 /h indicating robust QHE. In the lower inset, a 
more realistic magnetic field is chosen (45 T) to illustrate 
(i) the performance of the algorithm to probe low mag¬ 
netic fields and (ii) the possibility to destroy conductance 
quantization at higher energy for W = 1. For such dis¬ 
order, only the zero-energy Landau level fully develops, 
while all states become localized for E > E\. 

As another challenging example, we present the anal¬ 
ysis of a weakly correlated potential. We chose a glob¬ 
ally sublattice-symmetry breaking potential that is lo¬ 
cally disordered. In contrast to the above Anderson dis¬ 
order model, which shifts all on-site energies at random, 
the present potential is only locally present, i.e. only 
for a fraction of sites that are selected at random. The 
disorder is chosen such that it breaks globally the AB- 
sublattice symmetry. For the results shown in Fig. []] 
(main frame), we use a strength of Vab = O .270 with 
p = 2.5% of sites affected. Besides ordinary conductance 
quantization at the values ± 2 e 2 /h, ± 6 e 2 /h,..., an addi¬ 
tional plateau is induced at zero energy with zero Hall 
conductivity.— The plateau is clearly visible for fields be¬ 
tween 9 T and 45 T. We find the width of the zero-energy 
plateau equal for all calculated magnetic fields. The cor¬ 
responding energy value of the onset of the plateau at 
a = 0 is determined by the strength of the potential 


which is given as pV ab = O.OO 570 (indicated by dashed 
vertical line). The transition to higher plateaus (i.e. from 
± 2 e 2 //i to ± 6 e 2 //i) is rather influenced by the magnetic 
field-dependence of Landau level position. 

It is interesting to study the emergence of such plateau 
when gradually increasing the concentration of impuri¬ 
ties. Results are plotted for the case of 9T in Fig. [I] (in¬ 
set). The plateau-length scaling with p is evident from 
the figure. For the considered energy resolution in the 
inset (77 = 0.0004), a trace of plateau onset is still visible 
for the lowest percentages of the impurity potential down 
to 0.5% (corresponding to pV ab = O.OOI 70 or 2.7 meV). 


IV. CONCLUSION 

To conclude, an efficient real space algorithm to com¬ 
pute the Hall conductivity with the Kubo formula has 
been presented and validated on simple graphene-based 
systems. Such approach should become a useful compu¬ 
tational tool to simulate QHE in very large size complex 
disordered materials such as polycrystalline graphene^, 
graphene subjected to weak van der Waals interaction 
such as by a boron-nitride substrate^, or other types of 
disordered two-dimensional materials. It should also al¬ 
low to corroborate the formation of zero-energy plateaus 
of <7 xy driven by disorder-induced critical states, and dis¬ 
connected from degeneracy lifting of Landau levels^, an 
issue of genuine fundamental interest in the context of 
topological interpretation of the quantized conductance. 
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